Search on a Hypercubic Lattice through a Quantum Random Walk: II. d=2 
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We investigate the spatial search problem on the two-dimensional square lattice, using the Dirac 
evolution operator discretised according to the staggered lattice fermion formalism, d = 2 is the 
critical dimension for the spatial search problem, where infrared divergence of the evolution oper- 
ator leads to logarithmic factors in the scaling behaviour. As a result, the construction used in 
our accompanying article 6] provides an 0{y/N In N) algorithm, which is not optimal. The scaling 
behaviour can be improved to 0{\/N In N) by cleverly controlling the massless Dirac evolution oper- 
ator by an ancilla qubit, as proposed by Tulsi [j. We reinterpret the ancilla control as introduction 
of an effective mass at the marked vertex, and optimise the proportionality constants of the scaling 
behaviour of the algorithm by numerically tuning the parameters. 

PACS numbers: 03.67.Ac, 03.67.Lx 



I. TWO-DIMENSIONAL SPATIAL SEARCH 

The spatial search problem is to find a marked object 
from an unsorted database of size N spread over dis- 
tinct locations, with the restriction that one can proceed 
from any location to only its neighbours while inspecting 
the objects. Classical algorithms for unsorted database 
search are 0{N), but quantum algorithms working with 
a superposition of states do better. The quantum spatial 
search problem was first formulated on a single hypercube 
in Ref. 1], and its two-dimensional version has received 
particular attention due to its critical nature [H-Q . In our 
accompanying article @, we have investigated the spa- 
tial search problem for d > 2 hypercubic lattices using 
the massless Dirac evolution operator, and obtained close 
to the optimal scaling behaviour of Grover's algorithm 
. We have also developed a broad picture of how the 
dimension of the database influences the spatial search 
problem, in analogy with statistical mechanics of critical 
phenomena. In case of d < 2, the evolution operator aris- 
ing from the massless Dirac operator is infrared divergent 
(as J d'^k/k'^ in the continuum formulation). That slows 
down spatial search algorithms by logarithmic factors in 
d = 2 0-|5|, compared to the optimal scaling form, and 
our algorithm suffers the same fate. In this article, we 
modify our algorithm by introducing an effective mass in 
the Dirac evolution operator, and demonstrate how that 
overcomes the infrared hurdle and improves the scaling 
behaviour. 

Discrete time algorithms for the spatial search problem 
are based on random walks. In quantum random walks 
Q, the state amplitude distribution is unitarily evolved, 
such that amplitude at each vertex gets redistributed over 
itself and its neighbours at every time step. Quantum 
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superposition allows simultaneous exploration of multiple 
possibilities, and this technique has become part of a 
variety of algorithms (Ref. [9] provides an introductory 
overview). 

We study the specific case of search for a marked ver- 
tex on a square lattice with N = vertices. In our 
algorithmic strategy ^Oj] , the free Dirac Hamiltonian, 



Hi,, 



V -I- /3m , 



(1) 



diffuses the amplitude distribution around the lattice, 
while the potential attracting the amplitude distribution 
toward the marked vertex provides the binary oracle, 



R = e-'^^ - 
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Upon discretising the Hamiltonian according to the stag- 
gered lattice fermion formalism [Till, the anticommuting 
Dirac matrices are reduced to location dependent signs. 



n-l 



(3) 



To construct discrete and local time evolution, we need 
to exponentiate the terms in the Hamiltonian such that 
the resultant unitary operators are local. The mass and 
potential terms are single vertex local terms, and so can 
be exponentiated easily. To obtain local exponential of 
the kinetic term connecting neighbouring vertices, we 
separate the "odd" and "even" parts of the Hamiltonian 
on the bipartite square lattice, and exponentiate each 
block-diagonal Hermitian part separately [l^ : 

i?trcc = H, + H, + H,n : U = cxpi-iHr) . (4) 

The 4x4 blocks of the Hamiltonian, corresponding to 
hypercubes with coordinate labels {0, l}**^, are (in terms 
of tensor products of Pauli matrices) 
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and = —H^ when operating on a hypercube with ah 
coordinates flipped in sign. The factors appearing in the 
Trotter's formula for the discrete time evolution operator 
then become 



Iterate Qs times 



Uo{e) ^ Cl - isV2Ho(e) 

with s = sin(r/V2) and |cp + 



(6) 



— 1. We choose the 



unitary quantum walk operator to be 

W = UeUo = exp(-^(i/o + He)T) + 0{t^) . (7) 

Our initial state for the spatial search problem is the 
unbiased uniform superposition state, |s) = \x)/\/N. 
For m = 0, our algorithm alternates between the oracle 
and the walk operators, yielding the evolution 



(8) 



Here t2 is the number of oracle calls, and ti is the number 
of walk steps between the oracle calls. Both have to be 
optimised, depending on the size of the lattice, to find 
the quickest solution to the search problem. 

With iterative unitary operations, spatial search algo- 
rithms produce results periodic in time. Unlike Grover's 
algorithm, however, the maximum probability of being 
at the marked vertex, P, does not reach the value 1. 
Augmenting the algorithm by the amplitude amplifica- 
tion procedure 13J, the marked vertex can be found with 
probability 0(1), and the overall complexity of the algo- 
rithm is then characterised by the effective number of 
oracle calls tij^fP- 



II. IMPROVING THE ALGORITHM 

We have argued Q that spatial search in d dimensions 
obeys two lower bounds, 



t2 > max{dAf^/'',7r\/]V/4} , 



(9) 



following from distinct physical principles of special rel- 
ativity and unitarity, respectively. The bounds cross in 
the critical dimension c? = 2, where logarithmic correc- 
tions to scaling behaviour are expected in analogy with 
critical phenomena in statistical mechanics. Such loga- 
rithmic slow down factors have been observed in earlier 
works 043) and we want to suppress them as much as 
possible by suitably adjusting the evolution Hamiltonian. 

When TO = 0, the quantum random walk provides the 
fastest diffusion, minimising t^. But to = also makes 
the evolution operator infrared divergent in two dimen- 
sions, which spreads out the amplitude distribution in the 
Hilbert space too much and decreases P . Introduction of 
TO ^ would slow down the diffusion, but would also 
regulate the infrared divergence through k"^ ^ k'^ + m^. 
For small enough to, the diffusion speed (and hence 
may not change much, but substantial modification of the 
contribution of \k\ Sj rn modes can alter the behaviour of 







— < 


» — 






R 







— 1 


1 — 






w 



FIG. 1: Logic circuit diagram for Tulsi's controlled quantum 
spatial search algorithm. R and W are the binary oracle and 
the massless Dirac walk operator, respectively. We use the 
generalisation with W — > W*^ . 



P. In such a case, an optimal value of to can be obtained 
by trading off the increase in t2 against the increase in 
P. For a finite lattice, the lattice size acts as the infrared 
cutoff, and so the optimal value of to is expected to be a 
function of the database size N. 



A. Tulsi's Algorithm 

Tulsi constructed an algorithm possessing the previ- 
ously described properties, by controlling the evolution 
operators using an ancilla qubit [4i]. His scheme is illus- 
trated in Fig.l, where the ancilla operators are 



cos S sin S 
— sin S cos S 



Z = 



-1 
1 



(10) 



Is) to the 



and the algorithm evolves the initial state |1) 
target state |(5) |t) with |(5) = x]|l). 

For S — 0, Tulsi's algorithm reduces to spatial search 
using the massless Dirac walk operator, which finds 
the marked vertex with probability Pq = 0(1/ In TV) 
using Qo = e{VN In TV) oracle calls 0, i]. Tulsi 
showed that with cos (5 = Q{l/\/\nN), the algorithm in- 
creases the probability of finding the marked vertex to 
Ps = 0(1) without changing the scaling of oracle calls 
Qs = 0{VN In TV) [4]. More exphcitly, the algorithm 
largely confines the evolution of the quantum state to a 
two-dimensional subspace of the iV-dimensional Hilbert 
space, whereby 



Ps 



1 



4 cos (5 



1 + (B^ - l)cos^5 



(11) 



(12) 



Here B = Bq = O(vhiTV) is a second moment con- 
structed from the eigenspectrum of W, and characterises 
the infrared divergence of the problem. We have also 
included explicit factors of 2°* that are appropriate for 
spatial search with staggered fermions, where different 
vertices of an elementary hypercube correspond to dif- 
ferent degrees of freedom, and essentially only the degree 
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of freedom corresponding to the marked vertex partici- 
pates in the search process [6] . 

The optimal value of the ancilla control parameter 
is obtained by minimising the algorithmic complexity 
(where factors of 2'' cancel) : 

(cos ,5)opt = {B^ - 1)"'/' ~ II B , (13) 




B. Rephrasing of Tulsi's Algorithm 

Tulsi's algorithm is not presented in the Dirac operator 
language — rather the ancilla qubit is introduced in an ad 
hoc manner. Here we relate its ingredients to the well- 
known properties of the Dirac operator, in order to gain 
some physical insight in to its dynamical behaviour. 

Consider two species of Dirac particles: one with only 
the mass term in the Hamiltonian (completely at rest), 
and the other with only the kinetic term in the Hamil- 
tonian (fully relativistic) . Associating the species index 
with the ancilla value, we have 

4cc = |0)(0|®il,„ , 

4cc = \l){l\®{Ho + H,) . (15) 

Now, pick m to give the farthest unitary evolution, i.e. 

(The sign provided by /? does not matter in this case.) 
With these choices, the two species evolve independently 
with maximum evolution contrast. The total Hamilto- 
nian then yields the second half of the iteration in Tulsi's 
algorithm — schematically, 

\{) w ) 
= {ciW)Z ^ Z{ciW) . (17) 

The first half of the iteration in Tulsi's algorithm is a 
controlled oracle, conditioned on the state \S), i.e. 

(csR) ^ xliciR)Xs . (18) 

For the vertices x ^ without the potential, it is the 
identity operation. That lets the "|0)" species remain 
at rest while the "|1)" species diffuses at full speed, and 
there is no mixing between the two. A consequence is 
that the amplitudes |0)(X)|x ^ 0) do not change from their 
initial value zero. They neither mix with the amplitudes 
|1) |x 7^ 0) nor get any contribution from the amplitude 
|0) (8) 1^ = 0) by a walk step. 
The potential at a; = couples the two species, as per 

icsR)s=0 = ZX2S = ZX2S-n ■ (19) 



As a result, when the diffusing "|1)" species reaches the 
marked vertex, part of it gets converted to the "|0)" 
species and stops diffusing. At the next iteration, part 
of the "|0)" species gets reconverted to the "|1)" species 
and starts diffusing again. The net effect is that the con- 
versions reduce the number of W operations for walks 
when they pass through the marked vertex. The state 
|0) (El \x = 0) thus acts like a trap, and the concentration 
of the quantum state amplitude at the marked vertex in- 
creases Ps- Since the conversion between species pauses 
the walk, it can be interpreted as an effective mass, but 
this effective mass is unusual in the sense that it appears 
only at the marked vertex. 

The determination of the optimal value of S requires 
an analysis of the^eigenspectrum of the operator W. As 
shown by Tulsi [4i], when the evolution is largely con- 
fined to a two-dime nsion al subspace of the Hilbert space, 
(cos(S)opt = e(l/Vln7V) and 6 « tt/2. Note that for 
6 = tt/2, Ea.p9| implies that the two species decou- 
ple completely. So the mixing of species per iteration is 
tiny. Also, the target state is essentially the trap state, 
\S) ^ \x — 0) |0) ® \x — 0), which is reached from 
the initial state |1) ^ \s) with 8(1) probability by ac- 
cumulating the mixing of species over many iterations. 
The trapping of the quantum amplitude, resulting from 
an effective mass at a specific location, is a noteworthy 
feature — it can be an important ingredient in physical 
applications of the algorithm. 



III. SIMULATION RESULTS 

We carried out numerical simulations of our quantum 
spatial search algorithm, with a single marked vertex, 
and both with and without ancilla control. The choice 
to keep the walk matrices {Uo and Ue), as well as the 
ancilla operators (Xg and Z), real was convenient for 
numerical simulations. 

As in the case of our d > 2 simulations we first 
scanned for the best values of the parameters s and ti to 
optimise the algorithm. We again found that correlated 
s — ti pairs simultaneously maximise P and minimise the 
corresponding value of ^2- With increasing ti, P increases 
somewhat and t2 decreases slightly, but they are minor 
improvements. The major difference was observed in the 
dependence of the optimal parameters on the lattice size. 
Without ancilla control and for fixed ti, the optimal P 
decreases and the optimal s increases with increasing L. 
The variable 9 = \f2tx sin ^ s — t\T is somewhat larger 
then TT (the value found in case of d > 2 ^]), and increases 
with increasing L. This dependence on the lattice size is 
a consequence of the infrared divergence, whereby all the 
spatial modes do not contribute to the search process 
with equal strength. On the other hand, with ancilla 
control and for fixed t\, the optimal values of P and 
s show little dependence on L. The variable Q is still 
a bit larger than tt, but it is more or less a constant. 
These features indicate that the infrared divergence of 
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TABLE I: Fit parameters for peak probability without ancilla 



s 


tl 


L 


ai 


bi 


Error 


1 

n/2 


3 


512—16384 


2.607 


-12.76 


2.42x10"^ 



the spatial modes is brought under control. 

For concreteness in further analysis, we stuck to the 
parameter choice s = 1/V2 and ii = 3, which is close to 
the optimal choice, both with and without ancilla con- 
trol. It is not easy to theoretically estimate B, and hence 
determine the ideal proportionality constant between the 
control parameter cos 5 and 1 / Vln N. To figure that out 
numerically, we performed simulations with three differ- 
ent choices: cos S = y^l/\nN, i/i/bJV, i/S/lniV. 

Our results for the dependence of the peak probability 
on the lattice size are displayed in Figs. 2 and 3. As an- 
ticipated, P decreases logarithmically as 1/ In iV without 
ancilla control, and asymptotically approaches a constant 
with ancilla control. Even the behaviour of the sublead- 
ing correction changes from 1/ In L without ancilla con- 
trol to 1 /L with ancilla control, reconfirming that ancilla 
control indeed eliminates infrared divergence of P. The 
values of our fit parameters are presented in Tables I and 
II, where error refers to the root-mean-square (rms) de- 
viation of the data from the fit. In Table II, we have also 
included the related values of Bg following from Eq. ([TT|) . 

Next we display our results for the dependence of the 
number of oracle calls on the lattice size in Fig. 4. As 
expected, ^2 increases with decreasing cos S, as species 
conversions slow down the walk. But this slow down is 
only by a multiplicative factor, and vN\nN asymp- 
totically approaches a constant in all cases. Furthermore, 
the subleading correction is well parametrised as 
suggesting that ^2 is less affected by infrared divergence 
of the problem than P is. There is some oscillatory pat- 



TABLE II: Fit parameters for peak probability with ancilla 
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tl 


L 


cos 5 


ai 


bi 


Error 


Bs 








^1/lniV 


0.2243 


0.873 


2.50xl0~* 


1.056 


1 


3 


256—8192 


^4/ IniV 


0.1717 


2.536 


2.47x10"* 


1.207 








1/8/lniV 


0.1321 


1.429 


1.53x10"^ 


1.376 



tern in the data [14| . while the approach to the asymp- 
totic behaviour becomes smoother with decreasing cos 5. 
The values of our fit parameters are presented in Table 
III. Again error refers to the rms deviation of the data 
from the fit, and the related values of -B^^o following from 
Ea. ([TT|) are included. 

For a more elaborate comparison with Tulsi's analy- 
sis, we check whether our fit parameters obey Ea. (|lip 
parametrised in terms of the second moment B. With- 
out ancilla control, values of P and t2 should be related 
by (4ai)"^/^ ~ 802 /tt. That leads to the comparison 
0.31 O 0.36, which is quite reasonable. Stating it an- 
other way, our estimates of B vary from 0.31 ^log2 N to 
0.36v/iofeiV. Possible reasons for the discrepancy are: 
(a) contribution of subleading terms neglected in Tulsi's 
analysis (e.g. from the states outside the two-dimensional 
subspace used in the evolution), (b) the notorious diffi- 
culty in accurately extracting parameters from asymp- 
totically logarithmic fits. 

With ancilla control, the values of B^ in Tables II and 
III match well; control over infrared divergence definitely 
helps in extraction of the scaling behaviour. They give es- 
timates of B varying from 0.27A/iog^~/V to Q.Zl^J\og2 N. 
As per Eq. ([T3|l . therefore, the optimal choice for the an- 
cilla control parameter would be 

(cos(5)opt « 3.5^1/ log2A^ « ^8-5/ In A^ ■ (20) 
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FIG. 2: (Color online) Peak probability at the marked vertex 
as a function of database size without ancilla control. We used 
tl = 3, and the curve is the fit Plogj A'' = ai -|- (fei/logj N). 
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FIG. 3: (Color online) Peak probability at the marked vertex 
as a function of database size for different values of the ancilla 
control parameter cos 5. We used ti — 3, and the curves are 
the fits P = ai + (bi/L). 
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TABLE III: Fit parameters for search oracle calls for various TABLE IV: Fit parameters for the algorithmic complexity for 
values of the ancilla control parameter cos S various values of the ancilla control parameter cos S 



s 


tl 


L 


cos 5 


ffl2 


&2 


Error 


Bs 




s 


tl 


L 


cos 5 


0-3 


63 


Error 






512- 


-16384 


1 


0.1412 


2.755 


2.25x10" 


-3 






1 






y/1/lnN 


0.7336 


-13.06 


5.31x10"^ 


1 

V2 


3 


512- 


-8192 


^l/luAT 


0.3463 


-2.782 


2.41x10" 


3 


1.059 




3 


1024—8192 


1/4/ IniV 


0.4911 


-24.30 


4.03x10"^ 




512- 


-8192 


1/4/ In TV 


0.2030 


-8.977 


1.54x10" 


3 


1.242 










y/8/lnN 


0.4207 


31.24 


1.53x10"^ 






512- 


-8192 


^8/ \nN 


0.1562 


3.290 


3.81x10" 


3 


1.351 

















IV. SUMMARY 



Our numerical simulations had database sizes varying 
from N = 2^^ to 2^^. So in order to maintain cos 6 < 1, 
we had to r estrict t he proportionality constant between 
cos S and -y/l/ IniV. Nevertheless, our largest parameter 
choice, cos S — y/8/\nN, is close to optimal. 

Finally, we combine the results for P and t2, to look 
at the scaling behaviour of the algorithmic complexity 
t2/VP- Without ancilla control, the effective number 
of oracle calls scales as ^/NlnN, with a-ij ^/a{ being the 
proportionality constant. Our results with ancilla control 
are displayed in Fig. 5, together with the fit parameters 
in Table IV. The scaling of the effective number of or- 
acle calls is improved to -s/^VhiiV, with the proportion- 
ality constant 03 essentially the same as 02/^/0!. Due 
to the oscillatory pattern in the data and non-negligible 
subleading corrections, our results are inadequate to nu- 
merically optimise 03. On the other hand, our estimates 
of B and Eq. p^ give (a3)opt ~ 0.45. That is consistent 
with the value for our close to optimal parameter choice 
cos (5 — \J%/ hiN , and hence we infer our best algorith- 
mic complexity to be t2/\/P ~ 0.45 -^/iV logj N. 



For the spatial search problem, d = 2 is the critical di- 
mension where infrared divergences appear. Introduction 
of a mass term in the evolution operator can suppress 
infrared divergences, and we have reinterpreted Tulsi's 
ancilla control of the spatial search algorithm as intro- 
duction of an effective mass at the marked vertex. 

Our numerical results demonstrate how ancilla control 
improves the scaling behaviour of the spatial search al- 
gorithm in d = 2. In particular, they agree with Tulsi's 
predictions 0], and validate his analysis criterion that 
the evolution of the quantum state is largely confined to 
a two-dimensional subspace of the A^-dimensional Hilbert 
space. The change in scaling of P with ancilla control is 
a clear signal for suppression of the infrared divergence. 
Asymptotic behaviour of <2 does not change, however. 
It retains the Vln N factor beyond the lower bounds 
in Eq.®, indicating that some effect of the critical be- 
haviour survives. Generically, logarithmic factors cannot 
be fully eliminated in critical dimensions for interacting 
models of statistical mechanics. Still it is an open ques- 
tion, whether the logarithmic factors found in Tulsi's al- 
gorithm correspond to the minimal extra cost to be paid 
for the d = 2 spatial search problem, or whether they can 
be reduced further. 
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FIG. 4: (Color online) Number of oracle calls as a function 
of database size for different val ues of cos 5. We used ti — 3 
and the curves are the fits tij ^ N log2 A'' = 02 -f ihilV). 
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FIG. 5: (Color online) Effective number of oracle calls as a 
function of database size for different values of the ancilla 
control parameter cos 5. We used t\ = ?> and the curves are 
the fits ta/i/PAloiTlV = as + (63/i). 
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